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FOREWORD 

This report explains the digital computer program used to determine the natural Ire 
quenches, mode shapes, forces, and moments for propellant sloshing in tanks of 
arbitrary shape. The description of this analysis has been prepared under the re- 
quirement of Part 1 • Section A. 1. of Contract NASS-1119J and satisfies in puit the 
provisions of this requirement. The remaining documentation requirements of the 
subject contract are fulfilled in the three companion reports GD IA-DDE64-059. 

CD |A - DDE 64 -060. and GD lA - DDEO-1-061. 
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SECTION I 


INTRODUCTION 


The hydrodynamic forces and moments are derived (in Reference 1) for tanks pos- 
sessing a longitudinal axis of symmetry. These forces and moments are given in 
terms of coefficients which depend only on the tank geometry. 

This report explains the steps used to obtain these coefficients, given the tank geom- 
etry, and the procedures used in the program check-out. A description of the routines 
used in the program is included as well as instructions for use of the program. The 
output of the digital program gives the spring-mass parameters associated with the 
system . 


The digital routines were developed and programed by Roger Barnes. 
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SECTION 2 

SYNTHESIS OF EQUATIONS 

2. 1 HYDRODYNAMIC EQUATIONS . The equations which describe the sloshing of 
liquid propellant in a tank consisting of an arbitrary curve rotated about an axis of 
symmetry are derived in Reference 1 and are: 


) 


t, T (t) sin 9 * (r,z) (2. 1) 

^ n n 

n 1 

F ' -M a (2.2) 
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n n n 
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(2.4) 


where 


T ; is the surface wave height 

F 7 and F 7 are forces in the x , and x s directions (see Figure 2-1) 
3 2 3 2 

T 7 is the moment about the x ^ axis 

s t f r . z ) a r e e i ge n f u n c t i on s 
n 

and 


£ (t) is given bv the solution of (2.5). 
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£ (t> * Of ft) • -K b o (t) - K 
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Fa “ 

The K are nondimensiona! frequencies. 1. is the distance from the center of 

n 0' . ( 

gravity of the liquid *o the undisturbed true surface, and L is the distance from the 
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and 
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M is the total mass and V the total volume of the liquid. 
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( 2 . 10 ) 


where A is to be determined from the equation 
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A 
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A M 

" n n 
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(M n F(z)) cos /i (z 


■» 


(M F(z)) sin fi (z + d) 
n n 


2z 

— . -d z ' L (2. 11) 

\T 


To find an approximation for tj. truncate the series to include N - 1 terms. Evaluate 
(2.11) at N + 1 points (z.) on the boundary to obtain the following set of linear equations 
for the determination of the A n , i.e. 


A 
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1,2,..., N + 1 (2. 12) 


where r F(z) is the equation of the boundary AB. 


2.2 SYNTHE SIS F O K DIGITAL SOLUTION. The boundary value problem for the 

determination of the eigenfunctions. 4 , and eigenvalues. K , is 

n n 


•> ') 
r< £ , 4 4 :T<L 

n 1 n n n 


r r ■/. 
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(2. 13) 
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(on AB) 


(2. 14) 


4 K 
n_ n 

o7. 1. n 


(on BC) 


(2. 15) 


It has been shown (Reference 2) that the K / I. may be found by minimizing' the 


_K UA 

I, 


(-) r 


/|^(r, B )| 2 


(2. 16) 


To reduce the magnitude of the variables in integration, make the problem described 
by (2. M) through (2. 16) nondimensional by letting 

R r a 

Z z a 

Ur (Ha. /a) (K//.) 


a the distance from B to the z-axis. 

If the free surface of the liquid does not intersect the center line of the tank. A and 
0 coincide. In this situation, let the distance from the z-axis to point C be designated 
b v R , , and define c R a . 


2 R \{ 2 : 

■it n z 


(2. 17) 


(on AB) 


<2. IK) 


(on BC) 
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where the integrations are over nondimensional limits (see Figure 2-2). 



Figure 2-2. Nondimensional Geometry 

Following the pattern of Reference 2. let ; be expressed as a linear combination of 
the eigenfunctions for a shallow tank and eigenfunctions for a deep cylindrical tank, 
and use a Rayleigh-Ritz technique. Ret 

A — (2 . 

1 . 


and substitute 
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into (2 


2.20); differentiate with respect to c . and set the result equal to zero to obtain 


N* 

X! 1 mn mnj mn 
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(2.23) 


where 
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The set 


of functions : * (R.Z) are chosen, for reasons mentioned previously . as 


- * (R. Z) 

n 


2n~ 1 


R 


n = 12 


M • N* 


j (Z - 1. a) 




n = M ■*- 1 


2. 26) 


where 


j 1 o)» u 


(2.27) 


and the first root of (2.27) is called j M + ^ 


The b are calculated as 
mn 
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j e J (j e J ( f) - fJ j c) J j e) 
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where 


In 3 
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Substitute the value of ? * from (2.26) into (2.24) to obtain the following expressions 

for a . Thus, by Green's theorem 
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(see Reference 3. p. 239) 
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In a similar manner 
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Once the b and a are determined the eigenvalue problem can bi sohtd >u.ldin„ 
ran ran k 

and the cor responding e n . 

The k^ 1 mode is given by 


M-t 

R Ee k CMR Z) 
L - J n n 
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Its dimensional form is 


N* 

\ (r. z) = £ 



(r/a. z/a). 


n-1 

It is this value which must be used in the evaluation of (2. 6), (2. 7) and (2.8). 
are given by 


(2.37) 
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Finally consider the evaluation of 1^ • 
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SECTION 3 

ROUTINES OF THE COMPUTER PROGRAM 

The routines of the computer program are coded in the FORTRAN IV language, with 
the exception of parts of a SHARE simultaneous equation routine (SOLVE) and a 
SHARE eigenvalue routine (RWEG2F). 

3. 1 DRIVER ROUTINES . The boundary value problem requires solution of the 
eigenvalues of (2.23). Elements of the a mn and b matrices are computed for ten 
eigenfunctions, of which the first five are polynomials. The integration to evaluate 
the a elements is described in Section 3.2. The Bessel functions needed for the 
elements of both matrices are computed by a routine described in Section 3.3. The 
resultant eigenvalues, and their eigenvectors, are found by the Jacobi method using 
the routine RWEG2F. 


Routine RWEG2F solves an eigenvalue problem of the form: 


\ 


X 


X 


B 


( 3 . 1 ) 


where the B matrix must be positive definite 
die 


A 


not satisfy this requirement, though the a matrix did 


Now it was found that the 
Therefore 


b mn matrLx 
. letting 


mn 


and 


B 


mn 


(3. 1) is rewritten: 


b 

mn 


X 


A' 



(3.2) 


where 


A" (15.3) 

The \ ' eigenvalues are found and then inverted to give the desired eigenvalues of 
(2.22). Fach corresponding eigenvector is normalized by its largest element. 

For each mode of oscillation entered as input data, the nondimensional frequency (K fl ) 
is computed by (2. Mb). Force and moment coefficients (> n . b^. and h^) are computed 
using (2 . 29) through (2.42), where the h * factor of h^ employs a Bessel function, 
Section 2.2. and integration. Section 2,2. The above coefficients arc then combined 
as in (2.47) through (2.22) to yield the spring-mass parameters (m , * , K 
m , £ . and 1 ). The center of rotation is assumed to be at the bottom of the tank. 
Therefore. L, d in the expressions for *■ . i , and I . The I * term of I is iound 

by (2.9). The first two terms of I arc found by the integration described in 
Section 2.2. 
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Calculation of the third term of 1, , requires the solution of a set of simultaneous 
linear equations, yielding the coefficients. A„. Equation (2. 121 * evndunted a, 
ted points on the rigid tank wall below the undisturbed tree surface. 1 he total num . 
of points selected, n, must be such that. 


n it R 


(1) 


max 


40.0 


where R is the maximum radius of the rigid tank wall belt • 
max 

surface . 


(2) n 


50 


Restriction (1) is necessary to pn 


revent the modified Bessel function in the equation 


ms of the equations completely insigniii- 


from exceeding about 10 1(> . making other ter .... - 

cam Restriction (2) is set by program core storage limitations. In gene tul.p - 

selected by the program will lie end points of segments with equally spaced points u 
keen . point It. and point C (if < - (» on Figure 2-2. The resultant set of simultane- 
ous linear equations is -solved by the routine SOI.VK. Integration, agam as m 
Section 5.2. is then done on (2.4«) with the coefficients. A^. 
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3.2 INTEGRATION. Integration is performed using Gauss’s Quadrature Formula: 



f(u) du = (b-a) 


.1. 

f 

1/2 


1/2 


t (t) dt 



ft 


s <t t 

- i - > 


s 

n 



a) 


(3.4) 


where n is taken to be 16 and u - (b-a) t - a 


The sum in (3.4) is taken by adding the smallest term to the next largest and so on. 
to obtain the most accuracy when f(u) is increasing from a to b. 


Two integration routines were written: 


1. INTG1 - - for line integrals of the form 


B 

/ 


f </•) dz 


A 


The line integral is taken in a counterclockwise direction along each segment of 
the rigid tank wall from A to B. (See Figure 2-1.) 


B 


Z9 


z 3 


f f (z) dz = J f (z) dz J 


(z) dz 


U 

• j f (z) dz 


(3.5) 


Z 1 


ZO 


z n-l 


where 


^ ^ /.) ? ^ z n-\ are 7 ~ c oordinab of end points oi segments entered 

as input; and z n is the /.-coordinate at tl »e outer radius oi the undisturbed 
free surface. 

2. 1NTG2 — for line integrals of the form 


J f(r 


f(r) dr 


The line integral is taken in a eounterelockwise direction around the tluid ol the 
crosb section of the tank as in Figure 2-1. 


3-2 
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B 


C 


f (r) dr -j f (r) dr + J" 
A B 


f (r) dr 


(3.6) 


With r-coordinate limits of integration for each segment of the rigid tank wall 
from A to B and the undisturbed free surface, (3.6) is rewritten 



C 


f (r) dr = 


r 2 r 3 a 

f f (r) dr - J f (r) dr - . . . ~j f (r) dr 


ri 


T9 


r n-l 


Rn 


f ' * 


(r) dr 


(3.7) 


where 

r 1? r*2, r3, . . . , r n _i are the r-coordinate of end points of segments 
entered as input; a is the outer radius; and Rj is the inner radius of the 
undisturbed free surface. 


Equation (3.4) is then applied to each term of (3.5) and (3. 7). 


3.3 BESSEL FUNCTIONS . Two types of Bessel Functions are computed; the first 
kind of order one (Ji (x) ) and its derivative (Jj (x) ), and the modified of order one 
(Il(x) ) and its derivative (l{ (x) ). 


The BESSEL routine computes, in double precision, J^(x) and Jj(x) by the ascending 
series; 


l-9 k (ff W 

k! (LI)! 
k=l 


(3.8) 


k i x \2k 

<-D (“ ) 
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(k ! )“ 


k 1 
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tel 


~l) k (ff k 

2 k ! (k+1) 


(3.9) 


3-4 



G D | A-DDE64-062 


Each term of the above series is found by multiplying the previous term by the appro- 
priate factor: 


(-) 

(- 1 ) ' 2 1 

for J (x) and the second part of J (x) (3. 10) 

(k) (k~l) 1 1 


(- 1 ) 



for the first part of J ' (x) 


(3.11) 


Since 0 x ' 14.863588 in the program, the above technique is sufficiently accurate 
for the range of x. 


Routine f’M >PBS computes, in double precision, the l ^ (x) and t , ' (x) f>y polynomial 
approximations, with t x 3.75. 

For -3 . 75 x 3. 75: 

I (x) - x ' 1/2 *• 0. 87890534 t“ 

4 6 

0. 51498869 t • 0. 15084634 t 

8 10 

* 0. 02658733 t - 0.00301532 t 

12 - 

* 0.00032411 t ; (3.12) 

M*) 

i' (x) = I (x) (3.13) 

1 O X 

where 

•;> 

I (x) 1 ’ 3.5156229 r 

o 

*3.0899424 t* * 1.2067492 f G 

8 10 
*0.2659732 t 0.0360768 1 

12 


•0.0045813 t 


(3. 14) 
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For 3.75 ' x 


I (x) x" 1/2 c X 0. 39894228 - 0.03988024 t ' 


- 0. 00302018 

_•> 

1 

-3 

+ 0.00163801 t 


“4 

—5 

- 0.01031555 

t 

+ 0.02282967 t 

- 0.02895312 

-6 

t 

- 0.01787654 t 




- 0.00420059 

t 



(3. 15) 


I, (x) = I o W 


1 (X) 


(3. 16) 


where 


-1/2 x 

1 (x) - x e 

o 


0.39894228 * 0.01328592 t 


- 1 


-9 -3 -4 

*• 0.00225319 t “ - 0.00157 505 t • 0.00910281 t 

- 0.02057700 t ° * 0.02635537 t 6 - 0.01047033 t 


0.00392377 t 


-8 


(3. 17) 


3-0 
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SECTION 4 

UTILIZATION OF COMPUTER PROGRAM 


4. I DESCRIPTION OF PROGRAM INPUT , Three sets of numerical data are re- 
quired as input for a computer run. Sections 4.1.1 through 4. 1 . 3 describe this data, 
which is subject to limitations that are outlined in Section 4. 1.4. All numbers are 
entered on HO-column input forms. Each digit of a number or' a decimal point (if 
needed) occupies one column. Unless otherwise specified in the following sections, a 
number must be entered with a decimal point. The decimal point may be anywhere 
within the columns allocated for that number on the input form. 

A number that must not have a decimal point must be right-adjusted. That is. all 
digits must occupy the right-most columns of those allocated for the number: there 
may be no blank columns between the number and the last column, inclusive. 

Extremely large or small numbers may be entered as a coefficient times ten to a 
power: 


J i i l'A. • AL iL-iii j l 


Here a decimal point must be in the coefficient, and the last digit of the exponent must 
occupy the last column of those allocated for the number. 


All coordinates required as input must be from a system in which the r -coordinate of 
the centerline is zero, and the /-coordinate increases upward. 
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4. 1. 1 Problem Input . The first line of the input form contains that set of data re- 
quired to set up the problem for the computer. 

Columns Data 

i through 10 Enter the number of segments that describe the rigid tank 
wall. The number must be right-adjusted, without a dec- 
imal point. 

1 1 through 20 Enter the number of modes of oscillation desired. The 
number must be right-adjusted, without a decimal point. 

21 through 30 Enter the liquid density (in pounds per cubic foot). 

31 through 40 Enter the r- coordinate (in inches) of the beginning of the 

segments that describe the tank. 

41 through 50 Enter the z-coordinate (in inches) of the beginning of the 
segments that describe the tank. 

4. 1.2 Tank Ge ometry Input . Beginning with the second line of the input form, a line 
of the data below is required for each segment that describes the rigid tank wall. Seg- 
ments must Ik 1 ordered in a continuous, counterclockwise path around the tank. The 
number of these Lines of input must equal the number entered in columns 1 through 10 
of the first line (i.e. number of segments). 

Columns Data 

1 through 10 Enter the r-coordinate (in inches) of the end of the partic- 
ular segment, having proceeded along it in a counterclock- 
wise direction around the tank. 

11 through 20 Enter the z-coordinate (in inches) of the end of the particu- 
lar segment, having proceeded along it in a counterclock- 
wise direction around the tank. 

21 through 30 Ixnive blank for a straight line segment: otherwise, enter 
as follows, depending on the type of segment: 

Elliptical. Enter the r-coordinate (in inches) of the center 
of the ellipse. 

Circular. Enter the r-coordinate (in inches) of the center 
of the circle. 


4-2 
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Columns 

21 through 30 
(Contd) 

31 through 40 


41 through 50 


51 through 00 


(51 through 70 


Data 


Parabolic. Enter the r-coordinate (in inches) of the vertex 
of the parabola. 

Leave blank for a straight line segment: otherwise, enter 
as follows, depending on the type of segment: 

Elliptical. Enter the z~ coordinate (in inches) of the center 
of the ellipse. 

Circular. Enter the z-coordinate (in inches) of the center 
of the circle. 

Parabolic. Enter the z-coordinate (in inches) of the vertex 
of the parabola. 

Leave blank for a straight line segment; otherwise 
enter, depending on the type of the segment, as fol- 
low s: 

Elliptical. Enter the semimajor axis (in inches) of the 
ellipse. 

Circular. Enter the radius (in inches) of the circle. 

Parabolic. Enter the directrix (in inches) of the parabola. 

Leave blank for a parabolic or straight line segment: 
otherwise enter, depending on the type of the segment, 
as follow s: 

Elliptical. Enter the semiminor axis (in inches) of the 
ellipse. 

Circular. Enter the radius (in inches) of the circle. 


Leave blank for a circular or straight line segment; enter 
the amount of counterclockwise rotation from the norm 
(see below) of the ellipse if the- segment is elliptical, or of 
the parabola if the segment is parabolic. It may be left 
blank for an angle of zero degrees. 
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i.e. A toroid must be described by four segments of the same circle: 


z 



3 . The liquid level must not indicate a completely full tank. 

*1. There can be no more than ten modes of oscillation. 

7). Do not enter the centerline segment unless it is, in reality, a part ot the 
rigid tank wall . 

4.2 DESCRIPTION OF PROGRAM OCTPIT. Printed output on a computer run con- 
sists of the following: 

1. I'ank Geometry. A complete definition of all segments entered as input. 
And for each liquid level case: 

2. Liquid Level, The z -coordinate (in inches) of the undisturbed free surface 
in the tank. 

;L Mass of Liquid. The mass (in pounds) of the liquid contained in the tank, 

4. Center of Gravity. The z -coordinate (in inches) of the center of gravity of 
the liquid. 

-> 

5. I ' and its four terms (in lb-in. ). 

(L Eigenvalue Statistics. These, for each mode requested, are the eigenvalue 
and the eigenvector and its normalizing factor. 

7. For each mode requested, the coefficients (K . > , b and h ) ot the 

1 n n n n 

force and moment equations. 
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8. For each mode requested, the following parameters for the spring-mass 
analogy: 

m (in pounds) 
n 

i (in inches) 
n 

K ct (in lb in . ) 
n 3 

9. Also, the following spring-mass parameters: 

m (in pounds) 
o 

& (in inches) 

<> 

I (in lb-in. ~) 
o 

These parameters are a summation for all modes of oscillation entered as 
input data, and printed out in item 8 above. Therefore, if a one-mode 
analysis is to be used, only one mode should be entered as input data. 

4.3 SETUP FOR A COMPUTER RUN . A run on the computer requires a card deck 
containing system control cards, the program binary deck, and input data cards. 

The input data is punched on cards from the 80-column coding form described in 
Section 4.1. Arrangement of the deck will be as shown in Figure 4-1. 

The program is designed to be run under the IBSYS monitor system. For the IBM 
7094 computer, execution time estimates should be based on about 10 millihours per 
segment per liquid level case. 


4-0 




Figure 4-1. Card Deck Arrangement 
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SECTION 5 

CHECK-OUT OF THE DIGITAL PROGRAM 

The output of the computer program was checked against the experimental and (if 
possible) theoretical data available for tanks comprised of spheres, toroids, oblate 
spheroids, right circular cylinders, and concentric right circular cylinders. 

The volume and center of gravity of the fluid as well as the first two integrals in 
(2.9) may be easily checked out by hand-calculation for any tank configuration. The 
nondimensional frequencies, may be checked for any tank configuration that has 
analytical or experimental results available. The operations involved in finding 
these five quantities constitute the major portion of the program. The rest of the 
program simply manipulates these terms. Thus if these five quantities are calculated 
properly and the multiplications and additions required for the evaluation of the rest 
of the terms are checked out for a particular case, one is assured that the entire pro- 
gram is working properly. The only term which must be checked out in a different 
manner is the third term of (2.9). The boundary value problem for ip i (see Reference 
J, Equations 3.10 and 3.11) may be solved explicitly for specialized cases only. The 
shapes checked out here were for a right circular cylinder at various liquid levels. 

For a cylinder of radius 60 inches and liquid level 50 inches the calculated value was 
1 . 36 x 10®, the computer obtained 1.3779 x 10®. For a liquid level of 240 inches, 
note the comparison of I Q on page 5-6. Since ^ 1 is approximated using Bessel func- 
tions. the routine is not very accurate for elliptical or spherical tanks. However 
if the tank is "close to cylindrical" in shape, the output is good. That is, for the 
two Atlas tanks the routine calculates I properly. A rough check on this term for 
any tank shape is to compare the four terms comprising 1^ j • If the third term is 
of different magnitude than the other three, it is not usable. The approximation 
then used is that 1^ ' equals the sum of the first and fourth terms. 

The routines for a parabolic section were checked by using a cylindrical tank of 
radius 60 inches with a parabolic bottom. The parabola had a directrix of 5000 and 
thus approximated a straight line. The difference between the output for this shape 
and the output for a cylinder with a flat bottom were different only in the third signi- 
ficant digit. 
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Fi^iu'f l shn\\s tiir* (‘omnas'lson of the love/t three frequencies of a spherical tank 
(radius of oO inches) nil h a graph From Reference r>. In this ti^ure. whore R is the 
radius nf the sjiiicrc and h the liquid height, note ti*.al the lowest frequencies fall on 
the dark line while the next two sets of frequencies appear at times to ire closer to 
tne expert mental ratine than the theoretical results. 

The check for tin'* lowest frequencies for an obi ale spheroidal tank is given in Future 
4-2. 4 Tie lank used in this program had a semimajor axi - of 24 . 2d and a sc mi minor 

axis of 14. ! :i and was compared with me results in Relercnee 4. 

The second linear mode shown m Figure fL-2 is the t lur’d mode in Relercnee 4 since 
the second mode in Reference 4 is due to nonlinear ellects. 

The comparison o{' frequencies for a toroidal lank- with results from Reference (1, i; 
shown in Figure 1 ■ • Note Lhe good agreement lor the toe lowest modes. 

Hie frequencies for a rue tank, inner radius 10 iiu*hes and outer radius HO inehes. 
were- checked against the analvtieal solution lor a liquid height o{ >0 incites. (Fse 
must he made of the tables in Reference 7.) The results appear in the table below . 

AN A LYTIC A F COM!*! TFR 

K ^ n , P i i : > t ) . f 1 n I ( > 

K 2 , PPu 2 . UP 1 

K :L4> i. :i.4Pb 

• ) 

K l.Phu )Ahk> 

1 

hue proviso- A rntln output nuv also be win -eked out against the analytical solution? 

(Reference I p. :»- v ) or a right circular cxlindriea! tank. because ol a diilercnce 

m normalization in Lie compute* program arid the analytical work. the parameters 

b . a ad it ma\ net be directly compared. However the roinbinat ions vl>~ - b 
n n n ■ * . n i* u 

h t! . and > i i > j . 1 A v hit*!* appear in the spruu: mums parameters are independent of the 
normalization. ITus if the spr lay - mas s parameters cheek out (using results of 
Reference t. p. ;,- v i the *> !y % . h n a re being calculated properly. Il desired. Slow. - 

ever, these eoml i] na * i ons eeuid be checked out bv thomsrhr? , Flu* comparison oi 
parameters Is giver- in the follow mg table for a tank oi radius UP inehes and height < 

■ >pp era’s tilled fo a depth oi 4 m inches with a liquid ha\ mg a density f 4 1. 
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ANALYTICAL RESULTS COMPUTER OUTPUT 


K 


1 


K 


2 


K 


3 


y 

b 


1 

1 


i 

1 


m 

1 


V“3 


m 

o 

£ 

o 

M 


I 

o 


120 . 

3.682 

10.663 

17.073 

4 

3.729 x 10 

-4 

9.0965 x io 
-4 

13.259 x io 
174.91 

6 

0.30837 x 10 

9.4628 x io 3 
2.4059 x U) 
112.96 

6 

2.7143 x io 
0.14973 ■ 10 11 


120 . 

3.6823 

10.663 

17.073 

0.5969 

0.71901 

1.0486 

175.0 

6 

0.30841 x 10 
3 

9.4639 x 10 

6 

2.4059 > 10 
112. 95 

2.7143 x 10 

1 1 

0. 15081 x 10 
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a, in 


3.10 • 

6.43 Q 
13. 13 O 


1.0 


T (t- U> 

\ n n 
Spheroids ; 


50 
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APPENDIX 
SAMPLE CASE 
FOR 

AT LAS TANK 


A- 1 


650 PUNCH FORM - GENERAL PURPOSE 
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